Mechanical activation and non-local rheology in granular flows 
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Based on the idea that the motion of an assembly of rigid particles is mostly determined by 
the volume fraction, we propose a non-local constitutive relation for dense granular flows. This 
rheology is tested against discrete element simulations and calibrated using the spatial relaxation 
of the rheology observed under a homogeneous stress. We show that the creep flow obtained below 
yield conditions is quantitatively described by the same liquid rheology as the shear flow obtained 
above yield conditions. In the spirit of the experiment of Reddy et al. [l|, we investigate the influence 
of a shear band far away from the studied zone on its mechanical behaviour. The dominant features 
of the experiment are recovered both in the numerical simulations and in the analytical model. The 
later quantitatively fits the numerical data without any adjustable parameter. 



Granular materials belong to the class of amorphous 
athermal systems. Like foams [HQ, emulsions 0], sus- 
pensions 0-Q or metallic glasses Q, they exhibit a 
dynamical phase transition between static and flowing 
states. Analogously to phase transitions of thermody- 
namic systems, this jamming transition exhibits a diver- 
gence of correlation lengths 0-U 



revealing the pres- 
ence of non-local cooperative processes called dynamical 
heterogeneities |l2|. In order to describe the constitu- 
tive behavior of such systems, it is therefore natural to 
adopt the Ginzburg-Landauphenomenological approach 
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of phase transitions [Hi, ll3Hl5l | . The main issue is then to 
identify the relevant control and order parameters. Fol- 
lowing the now classical Liu-Nagel diagram for jamming 
transition (lH ]. it is usually assumed that the solid-liquid 
mechanical transition is controlled by the shear stress t 
HQ EH. which, once rescaled by its critical value, de- 
fines the yield parameter y. For a granular system under 
a fixed confining pressure P, one gets a Coulomb criterion 
y = fi^ 1 t/P, where fi c is the dynamical yield friction co- 
efficient. The order parameter is then the shear rate j 
or any combination of 7 and the stress components, like 
the fluidity 7/r. However, different experiments have 
challenged this Coulomb picture, which assumes a tran- 
sition between solid (7 = 0) and liquid (I7I > 0) states at 
y = 1. (i) In the inclined plane geometry, thin granular 
layers flow anomalously [17( and stop at a yield param- 
eter y > 1 [TH, EH- (ii) A creeping flow is commonly 
observed in regions which are expected to be jammed, 
since y < 1. [191421 (iii) A solid plunged in grains and 
submitted to a force lower than the yield threshold starts 
moving as soon as a shear band is created far away from 
the solid [H, [22|]. Moreover, from the mechanical point 
of view, the choice of y as a control parameter suffers 
from a major flaw: the stress is not amongst the state 
variables. Indeed, following Newton's law, the state of 
the system is determined by the particles positions and 
velocities, from which forces are deduced. To say it dif- 
ferently, a dynamical phase transition model based on 
the assumption that the stress field is known cannot be 
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FIG. 1: (Color online) (a) Relation between the effective fric- 
tion coefficient fi and the volume fraction <f>. (b) Schematic 
of the numerical set-up used to calibrate spatial relaxation 
to the bulk constitutive relation, (c) Typical profiles of the 
yield parameter y obtain numerically below (thick black line, 
J4| < 1) and above (thin red line, |34| > 1) yield conditions, 
(d) Corresponding velocity profiles. The dotted lines show 
the predictions of a local rheology. Pb is the pressure in SB. 

introduced in continuum dynamical equations to extend 
Navier-Stokes to yield stress fluids. To do so, one needs 
a formulation where stresses are deduced from true state 
variables: volume fraction, strain, strain rate, coordina- 
tion number, fabric tensor, etc. The aim of this letter 
is to propose a phenomenology for dense granular flows 
and to test it against discrete element simulations. 

Before proceeding to the derivation of the non-local 
rheology, we wish to emphasize that Ginzburg-Landau 
approach does not capture much (if anything) from the 
dynamical mechanisms at work, which (unfortunately) 
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constitute a separate issue. The stress-based approach 
implicitly assumes that a flow can be decomposed into 
a quasi-static sequence of plastic events. In this picture, 
localized plastic rearrangements get triggered when the 
local stress state crosses a microscopic local yield crite- 
rion ll|, |23|, |24j , namely, a local Coulomb criterion in 
the case of granular matter. These localized events in- 
duce long-range fluctuations of the stress over the sys- 
tem, which can "activate" a cascade of new plastic events. 
The mechanical fluctuations in amorphous athermal sys- 
tems would then play the role of temperature in ther- 
mal systems: rearrangements would be mechanically ac- 
tivated, suggesting a picture analogous to Eyring's tran- 
sition state theory for the viscosity of liquids [H, [23| . The 
shear rate 7 then gives the rate at which quasi-static re- 
alizations of the contact force network are generated. 

This kinetic elastoplastic picture is challenged by the 
geometrical picture, which has been developed for dis- 
ordered packings of particles that interact through re- 
pulsive contact forces. In this alternative approach, me- 
chanical properties are related to the contact network 
geometry and to the distance to isostaticity [j| [25 - 27 [ . 
It has been recently shown that the statistical properties 
of grain trajectories and in particular their cooperative 
non-affinc motions are mostly controlled by the volume 
fraction (j), whatever the nature of the dissipative mecha- 
nisms [28[ . To move by a distance as the crow flies equal 
to its diameter d, a grain makes a random-like motion 
whose average length diverges as ~ d((j> c — (j))" 1 , lead- 
ingto a stress tensor diverging as /(</>) ~ {<f) c — t/>)~ 2 
[fiTS HS|- Provided the flow is homogeneous, the stress 
ratio t/P is related to the volume fraction <p by a function 
p(4>) (Fig.[T^,) which does not depend on the interparticlc 
friction coefficient nor on the flow regime (overdamped or 
inertial) Q. In this geometrical picture, the "local plastic 
events" arc replaced by " non-local cooperative motions" ; 
the ability to flow is governed by steric effects and con- 
trolled by the volume fraction in the neighborhood of the 
particle. 

Following the geometrical picture, the lowest order cor- 
rection to a local rheology [ill [2t| involves the parameter 
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when the linearity and the isotropy of the response to a 
variation of 4> is assumed. The rheology of a dense inertial 
granular flow can then be written under the form: 
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where jij = ^ff 1 + j^jr^ is the rate of strain tensor and 
p the density, x is t ne non-local correction, which ver- 
ifies x( K ) — Xo K + 0(k 2 ) (Fig. (2b). K is positive when 
the point considered is surrounded by a lower volume- 
fraction: the shear stress therefore decreases. For a sim- 
ple shear, the volume fraction can be eliminated between 




FIG. 2: (Color online) (a) Relaxation length I below (o) and 
above (□) yield conditions. The solid line is a phenomenolog- 
ical fit to the data that diverges as 34 — 1| _1/ ' 2 . (b) Function 
x(k) measured for 34 < 1. The red solid line shows the slope 
Xo — 8 deduced from the best fit of the data for 34 > 1. 

the shear stress r = o~ xz and the pressure P — —a zz . 
The non-local rheology can then be expressed as: \y\ = 
(1 — x( K )) m(0)/Mc- However, formaly as //(</>) and /(</>) 
are linear [30j , the very same non-local rheology can be 
expressed in terms of the inertial number / = \j\d/ y/P/ p 
rather than 0, to match the kinetic elastoplastic picture. 
/ compares the rate at which plastic event occurs ~ I7I 
to the inertial time-scale ~ dj y/P/p over which these 
events take place. One gets: 
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In both formulations, we consider that the whole granular 
assembly flows, whatever the yield parameter 34 as soon 
there is a flow somewhere. The very same rheology then 
applies to fast flows, above the yield criterion y > 1, and 
to evanescent flows, below the yield criterion y < 1. 

We have calibrated and tested this non-local rheology 
using molecular dynamics simulations (see [3(| for de- 
tails). The system is two-dimensional and slightly poly- 
disperse. Contacts are frictional. The shear cell is com- 
posed by two rough walls, with periodic boundary condi- 
tions along the direction x parallel to the walls (Fig.[T}3). 
The position of the walls is controlled to ensure a con- 
stant normal stress and a constant velocity. We consider 
the limit of rigid grains for which results do not depend 
any longer on the microscopic spring constant nor on the 
damping time. The originality of the numerical set-up is 
the presence of bulk forces applied to the grains, which 
depend of their positions (Fig. [TJd). The system behaves 
as if there was gravity field varying in space. The bulk 
forces can be either along z, to vary P, or along x, to 
vary r. We are therefore able to impose the field of y, 
once the system as reached a steady state (Fig. []}:). 

We first impose homogeneous stress conditions in the 
center of the shear cell (noted SB for shear band in 
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FIG. 3: (Color online) a) Schematic of the numerical set-up 
used to investigate non-local effects, (b) Typical profile of the 
yield parameter y. (c) Velocity profile, (d) Velocity gradient 
profile. Pb is the pressure in both SB and SSR. 

Fig. Hb). The later is surrounded by two buffer zones 
hereafter called the master shear zones (MSZ) and in 
which y is gradually varied using vertical forces (Fig. [lb). 
One observes that the velocity profile in the central 
zone is not linear but is well fitted by the form jbZ + 
C sinh(z /£) as predicted by the linearization of Eq. ((5J) 
around the bulk inertial number (Fig. QJ1). Above yield- 
ing conditions, the predicted relaxation length £ reads 
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for \y b \ > 1 



V \y»\ - 1 

while below yielding conditions, one gets 7b = and 

for \y b \ < 1 
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Therefore, the relaxation length £ should diverge sym- 
metrically on both sides of the yielding condition |34| = 
I, as - ||;y 6 | - l^ 1 / 2 . Fig. Hi shows that this 

relation is nicely verified by the data. Let us emphasize 
again that both sides of the divergence correspond to the 
liquid granular state and are predicted by the very same 
non-local correction xi K )- F° r |34| < 1 the grains in the 
SB only flow due to the influence of the MSZ, which are 
above yielding conditions. Importantly, the existence of 
a linear response of the system around the critical state 
{1 = and <fi = cf) c ) imposes to build k by rescaling V 2 <fi 
by</> c -0(or V 2 /by/, Eq.E}. 

We now carry out a numerical experiment similar in 
principle to that performed by Reddy et al [l| in a Cou- 



ette cell. The conceptual idea is to measure the rheol- 
ogy in the shear band, which is below yield conditions, 
under the influence of the master shear zone (MSZ) lo- 
cated close to the rigid boundaries. To implement such 
a slaved secondary rheometer (SSR), the shear stress is 
changed in a small region of width 2H S around the cen- 
ter of the shear cell, by applying horizontal forces to the 
grains (Fig. [3]). It is characterized by the yield parameter 
y s and by the central shear rate 7 S . The flow in the SSR 
is slaved to that forced in the MSZ, which is character- 
ized by a shear rate j m . The SSR is separated from the 
MSZ by a region of width Hb characterized by a yield 
parameter < 34 < 1- So, there is a single observable, 
j s and five parameters: j m , 34, Hb, 34 and H s . 

Experimentally, the SSR used in [lj is a small im- 
mersed rod submitted to an external force, which, once 
rescaled by its critical value, plays the role of 34 • The 
rod creep velocity is the analogous of j s . The shear rate 
7 m of the MSZ is controlled by the rotation velocity of 
the inner cylinder. Note that contrarily to the numeri- 
cal set-up, the radial profile of the yield parameter y is 
inhomogencous in the Couctte cell. The three key obser- 
vations reported in [l[ are recovered here, (i) The shear 
rate j s in the SSR is proportional to the shear rate 7 TO 
in the MSZ (Fig. 2k). (ii) j s (roughly) decreases expo- 
nentially with the distance to the yield conditions in the 
SSR, measured by 1 — |34|- (Fig. Ek) (iii) j a decreases 
exponentially with Hb (Fig. [4)}). In [lj, the property (i) 
has been interpreted as the signature of a quasi static flow 
and (ii) as a Boltzmann-like factor suggesting an Eyring- 
like mechanical activation. This Boltzmann-like factor 
was itself related in [23| to the exponential tail of the 
contact force distribution, using the kinetic clastoplastic 
approach. 

Can the three properties shared by the experiment 
and the numerics be recovered in the non-local formal- 
ism derived here? Using again the linearization of Eq. [2] 
around the critical state, the solution in the SSR takes 
the form 7 = j s cosh [z/£(34)] while in the SB, it reads 
7 = 7+ cxp [z/l{y b )\ +4_ exp[-z/£(y b )}. At the in- 
terface z = H s of the secondary rheometer, the model 
predicts that \^\ is continuous even when 7 changes sign. 
This property is remarkably verified in the numerics 
(Fig. [3ji). Matching |4| at z — H s , one gets 
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in the limit Hb S> ^(34)- This expression is reported in 
Fig. |4] without any adjustable parameter. The agreement 
with the numerical points is almost perfect when the SSR 
is sheared in the direction opposite to the MSZ (34 < 0). 

The non-local model predicts the proportionality of j s 
and 7 m (property i) as a trivial consequence of the lin- 
earization of the rheological equation. It explicitly pre- 
dicts that the influence of the distance Hb can be factor- 
ized and is exponential (property iii), due to the spatial 
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FIG. 4: (Color online) Ratio of the shear rate in the secondary rheometer 7 S to that in the master shear band j m . (a) Linearity 
of the response, (b) Data obtained for different distances Hb collapse on a master curve, once rescaled by the factor giving the 
influence of the distance: 7a/ = 7 m exp (— Ht,/£(yb)). (c) Influence of the size H s of the secondary rheometer. The solid lines 
are the predictions of the non-local rheology obtained without any adjustable parameter. 



relaxation of 7 in the zone separating the SSR from the 
MSZ. The fast exponential-like decay of % with 1 — \y s \ 
(property ii) also results from the spatial relaxation of 
the shear rate, but this time inside the SSR. 

Our results demonstrate that the granular material 
constitute a single liquid phase across yielding conditions 
(y = 1), which continuously connects the fast shear flow 
(MSZ) to the evanescent flow in the secondary rheome- 
ter (SSR), even when the direction of the flow is reversed 
(continuity of I7Q. The three properties shared by the ex- 
periment, by the numerics, and by the analytics turn out 
to be simple consequences of the linearization around the 
critical state. Exponential factors are therefore neither 
associated to Boltzmann factors nor to the tails of the 
contact force probability distribution function, as previ- 
ously assumed. Finally, Fig. [4ji shows a blind prediction 
that can be easily tested experimentally, although not 
reported in [l|. One observes that the size of the SSR 
(~ the rod) has a strong influence on 7 S (~ the creep 
velocity): the wider the SSR, the faster the decay of 7^ 
with the distance to yield conditions. The non-local rhe- 
ology proposed here must now be included into Navier- 
Stokes solvers to test its predictions in more demanding 
geometries like the split-bottom setup or the channelized 
avalanches, as well as for time-dependent flows. 
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FIG. 5: a) Schematics of the simulation set-up. Walls are 
depicted by dark purple circles, (b) Same, but now also with 
the secondary rheometer illustrated. 

SUPPLEMENTAL MATERIAL 
Numerical model and viscoelastic parameters 

We consider a two-dimensional system constituted 
of ~ 2 • 10 3 circular particles of mass and diameter 
di, with a ±20 % polydispersity (Fig. The shear 
cell is composed of two rough walls distant by H and 
each moving along the ^-direction at different constant 
velocities, with the difference of the velocities equal to 
u% . The position of the walls is controlled to ensure 
a constant normal stress P w at the walls (whereby H 
will change/fluctuate during the simulations). Periodic 
boundary conditions are applied along x. The particle 
and wall dynamics is integrated using the Verlet algo- 
rithm. The particles are submitted to contact forces 
modeled as viscoelastic forces along the normal direction 
and as a Coulomb friction along the tangential direction 

M. 

The normal spring constant k„ is chosen sufficiently 
large (i.e. k n /P > I0 3 ) to reach the rigid asymptotic 
regime. The Coulomb friction coefficient is chosen equal 
to fi p = 0.4 and damping parameters are chosen such to 
yield a restitution coefficient e ~ 0.9. 
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FIG. 6: The volume fraction <j> as a function of the inertial 
number /. Numerical data points are given by triangles and 
Eq.[6]by the solid line. 



Where the origin of the z-axis is taken in the center of 
the cell, §t = (H — Hsb) an d &b = (Hsb — H), a ~ d 
and f z is the magnitude of the forcing. 

We furthermore introduced a second (stress controlled) 
rheometer in our simulation set-up. This was done by 
applying two horizontal anti-directional stresses within 
the zone below yield conditions separated a distance 2H S 
apart. The second rheometers imposed forces were given 
by distributing a stress over a horizontal layer of grains on 
respectively sides of the entrance of the second rheometer 
(as illustrated in Fig. 1). 



Analytical solution 



The local rheology of a granular flow is given by, 

<j>(l)=<j> c -al and M (0) = fx c (l + p(<j> c - <j>)) . (6) 

For the specific system with the values <fi c = 0.8153, [i c = 
0.2665, a = 0.3112, and p = 15.6127 (see Fig. 2 in the 
SI and Fig. la in the Letter). 

The nonlocal rheology in terms of the inertial number 
I (analogous derivation can be written in terms of the 
volume fraction (f>) is given by (Eq. 3 in the Letter): 



External forces and the second rheometer 

In order to create a zone below the yield condition 
(i.e. \y\ < 1), two anti-directional vertical gravity-like 
forces were applied on the grains with gaussian spatial 
distributions (on particle labelled i): 
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(7) 



where \ 1S the non-local correction, which to a first ap- 
proximation reads x( K ) — Xo k + 0( k2 )- We first consider 
the situation above yielding conditions |^| > 1. Lineariz- 
ing the equations around the zero order 1 3^ | = 1 + PIo, 
where I is the local approximation, one gets: 



x'o\y\ 



d 2 V 2 h = 0. 
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FIG. 7: Red solid line: The full nonlocal correction x as a 
function of k according to Eq. [10] Black dashed line: Linear 
approximation of Eq. 6 (i.e. \ ~ Xo K an d with Xo — 8). 



Once again one finds that I\ is given on the form ~ 
cosh(z/f), but now with I = , d . 

For small k (i.e. close to the yield condition \y\ = 1) 
the length t is found to be well approximated by ~ 
(iyXo \\y^\ ~ 1| _1//2 - For larger k and below the yield 
condition we, however, found that the non-local function 
X took the shape (see Fig. 3): 



y/ (I — na) 2 + nb(na — 1) 



an 



(10) 



with a = —15.95 and b = 16.3. This gives rise to an 
asymmetry in I around \$\ = 1 when sufficiently far away 
from the same point (as seen in Fig. 2 in the Letter). 



The solution Ii (i.e. the non-local correction to the in- 
ertial number) takes the form ~ cosh(z/£) with I = 




Below the yielding condition |^| < 1 we instead linearize 
the equations around the critical state Iq = and find: 

\y\=i- x ( K ) or x -\i-\y\) = d 2 ^±. (9) 
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